{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2b04f3113ae44cc68df17ba7e11b7f8f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Plot(antialias=3, axes=['x', 'y', 'z'], background_color=16777215, camera=[4.5, 4.5, 4.5, 0.0, 0.0, 0.0, 1.0, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"The Lotka-Volterra model where:\n",
    "x is the number of preys\n",
    "y is the number of predators\n",
    "\"\"\"\n",
    "#Credits:\n",
    "#http://visual.icse.us.edu.pl/NPB/notebooks/Lotka_Volterra_with_SAGE.html\n",
    "#as implemented in K3D_Animations/Lotka-Volterra.ipynb\n",
    "#https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "def rhs(y0, t, a):\n",
    "    x, y = y0[0], y0[1]\n",
    "    return [x-x*y, a*(x*y-y)]\n",
    "\n",
    "a_1 = 1.2\n",
    "x0_1, x0_2, x0_3 = 2.0, 1.2, 1.0\n",
    "y0_1, y0_2, y0_3 = 4.2, 3.7, 2.4\n",
    "\n",
    "T = np.arange(0, 8, 0.02)\n",
    "sol1 = odeint(rhs, [x0_1, y0_1], T, args=(a_1,))\n",
    "sol2 = odeint(rhs, [x0_2, y0_2], T, args=(a_1,))\n",
    "sol3 = odeint(rhs, [x0_3, y0_3], T, args=(a_1,))\n",
    "\n",
    "limx = np.linspace(np.min(sol1[:,0]), np.max(sol1[:,0]), 20)\n",
    "limy = np.linspace(np.min(sol1[:,1]), np.max(sol1[:,1]), 20)\n",
    "vx, vy = np.meshgrid(limx, limy)\n",
    "vx, vy = np.ravel(vx), np.ravel(vy)\n",
    "vec = rhs([vx, vy], t=0.01, a=a_1)\n",
    "\n",
    "origins = np.stack([np.zeros(np.shape(vx)), vx, vy]).T\n",
    "vectors = np.stack([np.zeros(np.shape(vec[0])), vec[0], vec[1]]).T\n",
    "vectors /= np.stack([np.linalg.norm(vectors, axis=1)]).T * 5\n",
    "\n",
    "curve_points1 = np.vstack([np.zeros(sol1[:,0].shape), sol1[:,0], sol1[:,1]]).T\n",
    "curve_points2 = np.vstack([np.zeros(sol2[:,0].shape), sol2[:,0], sol2[:,1]]).T\n",
    "curve_points3 = np.vstack([np.zeros(sol3[:,0].shape), sol3[:,0], sol3[:,1]]).T\n",
    "\n",
    "########################\n",
    "from vtkplotter import *\n",
    "\n",
    "Arrows(origins, origins+vectors)\n",
    "\n",
    "Line(curve_points1, c='o', lw=8)\n",
    "Line(np.vstack([T, sol1[:,0], sol1[:,1]]).T, c='o')\n",
    "\n",
    "Line(curve_points2, c='g', lw=8)\n",
    "Line(np.vstack([T, sol2[:,0], sol2[:,1]]).T, c='g')\n",
    "\n",
    "Line(curve_points3, c='b', lw=8)\n",
    "Line(np.vstack([T, sol3[:,0], sol3[:,1]]).T, c='b')\n",
    "\n",
    "show(collection()) # all sofar created objects\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
